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Abstract 

Background: The biological basis for sex differences in brain function and disease susceptibility is poorly understood. 
Examining the role of gonadal hormones in brain sexual differentiation may provide important information about sex 
differences in neural health and development. Permanent masculinization of brain structure, function, and disease is 
induced by testosterone prenatally in males, but the possible mediation of these effects by long-term changes in the 
epigenome is poorly understood. 

Methods: We investigated the organizational effects of testosterone on the DNA methylome and transcriptome in two 
sexually dimorphic forebrain regions — the bed nucleus of the stria terminalis/preoptic area and the striatum. To study the 
contribution of testosterone to both the establishment and persistence of sex differences in DNA methylation, we 
performed genome-wide surveys in male, female, and female mice given testosterone on the day of birth. Methylation 
was assessed during the perinatal window for testosterone's organizational effects and in adulthood. 

Results: The short-term effect of testosterone exposure was relatively modest. However, in adult animals the number of 
genes whose methylation was altered had increased by 20-fold. Furthermore, we found that in adulthood, methylation at 
a substantial number of sexually dimorphic CpG sites was masculinized in response to neonatal testosterone exposure. 
Consistent with this, testosterone's effect on gene expression in the striatum was more apparent in adulthood. 

Conclusion: Taken together, our data imply that the organizational effects of testosterone on the brain methylome and 
transcriptome are dramatic and late-emerging. Our findings offer important insights into the long-term molecular effects 
of early-life hormonal exposure. 
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Background 

The biological basis for sex differences in brain function 
and disease susceptibility is poorly understood. Numer- 
ous neurological disorders (e.g., autism, schizophrenia, 
Parkinson's disease) show sexual dimorphism in preva- 
lence [1-3], and sex-specific biological factors are likely 
to be major contributors. Sex steroid hormones such as 
testosterone play a major role in sexually dimorphic 
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brain development [4,5]. Exposure of neural tissue to 
testosterone and estradiol, its aromatized form, during 
the perinatal window ('the sensitive period'), leads to 
long lasting and irreversible masculinization of the brain 
[6]. The fundamental molecular mechanisms underlying 
the hormonal regulation of brain sexual differentiation 
remain understudied. 

Emerging evidence implicates DNA methylation as an 
important player in a variety of critical nervous system 
functions [7-9]. 5-Methylcytosine (5-mC) marks at CpG 
islands in gene promoters are known to affect gene tran- 
scription [10], modulate X inactivation and imprinting, 
and regulate heterochromatin [11]. There is also increas- 
ing evidence of the importance of non-CpG methylation 



© 2014 Ghahramani et al.; licensee BioMed Central Ltd. This is an Open Access article distributed under the terms of the 
Creative Commons Attribution License (http://creativecommons.Org/licenses/by/4.0), which permits unrestricted use, 
distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public 
Domain Dedication waiver (http://creativecommons.Org/publicdomain/zero/1.0/) applies to the data made available in this 
article, unless otherwise stated. 



Ghahramani et al. Biology of Sex Differences 2014, 5:8 
http://www.bsd-journal.eom/content/5/1/8 



Page 2 of 18 



in neural development [12,13]. Recent studies have iden- 
tified CpGs that show sex-specific methylation that can 
be modified by sex steroids during the sensitive period. 
Estradiol can alter the methylation status of CpG sites 
on the promoters of the estrogen and progesterone re- 
ceptor genes [14,15]. However, the methylation status of 
only a limited number of CpG sites at candidate genes 
has been studied. A larger-scale study of the methylome 
could elucidate the role of epigenetic modifications in 
hormone-induced brain sexual differentiation. 

We hypothesized that long-term effects of gonadal 
steroid hormones on brain development involve epigen- 
etic modifications. We investigated the scope and prop- 
erties of testosterone's organizational effects on neural 
DNA methylation and gene expression by comparing 
the methylome and transcriptome in male (XY), female 
(XX), and female mice that received testosterone on the 
day of birth (XX + T). Genome-wide methylation and 
expression profiling were carried out for two sexually di- 
morphic brain regions: the striatum and the combined 
bed nucleus of the stria terminalis and preoptic area 
(BNST/POA). These regions were chosen as both are re- 
sponsive to gonadal hormones and show strong sexual 
dimorphisms. The BNST/POA is subject to long-lasting 
irreversible neuroanatomical changes due to perinatal 
testosterone exposure and has been implicated in the 
regulation of male copulatory behavior, gonadotropin re- 
lease, and stress modulation [16]. The striatum is in- 
volved in dopaminergic function and reward and shows 
several key sex differences, many of which are caused by 
gonadal hormones. Numerous aspects of dopamine me- 
tabolism are influenced by estrogen [17-20]. Further- 
more, we have previously shown that Sry (the Y-linked 
male sex determination gene) enhanced striatum dopa- 
mine release and regulated sensorimotor functions of 
dopaminergic neurons [21]. To examine the long-term 
molecular effects of organizational testosterone, we ex- 
amined two different time points: postnatal day (PN) 4, 
which is during the sensitive period, and adulthood 
(PN60). 

We show that neonatal testosterone treatment of fe- 
males induces a shift in the methylome from a female- 
typical to a more male-typical pattern. Contrary to our 
expectations, the shift toward the male pattern is only 
observed during adulthood. Organizational testosterone 
also affects the CpG methylation status of many more 
genes at PN60 than at PN4. Consistent with the CpG 
methylation analysis, testosterone's effects on gene ex- 
pression in the striatum were more apparent in adult- 
hood. Our data demonstrate that the organizational 
effects of testosterone on the brain methylome and tran- 
scriptome are dramatic and late-emerging. They suggest 
an important role for CpG methylation in brain sexual 
organization. 



Methods 

Animals and neonatal injections 

This study was approved by the University of California, 
Los Angeles (UCLA) Committee on Animal Research 
and was performed in accordance with the recommen- 
dations in the Guide for the Care and Use of Laboratory 
Animals of the National Institutes of Health. C57BL/6 J 
female and male mice were purchased from Jackson La- 
boratories (Bar Harbor, ME, USA) and housed at the 
UCLA Animal Care Facility. Animals were maintained 
at 20°C with a 12-h light/ 12-h dark cycle, provided ad 
libitum with food and water, and allowed to acclimate 
for 1 week before initiation of experiments. Female mice 
were mated and once pregnant, cages were checked daily 
for pups. 

On the day of birth, PNO, male pups (referred to as 
XY throughout this paper) were treated subcutaneously 
with 15 [A of peanut oil vehicle; female pups were either 
treated subcutaneously with 15 \A oil (XX) or with 
100 \ig testosterone propionate ( Sigma- Aldrich, St. 
Louis, MO, USA) in 15 [A oil (XX + T). Half of the mice 
in each group were sacrificed at PN4 during the peri- 
natal sensitive period for sexual differentiation and just 
prior to the sexually dimorphic cell death known to be 
important for differentiation of the POA and BNST 
[22,23]. The remaining mice in each group were sacri- 
ficed in early adulthood, on PN60. In order to eliminate 
group differences in pubertal or adult gonadal hormone 
exposure, mice euthanized on PN60 were gonadecto- 
mized at PN21 (i.e., prior to puberty) and implanted 
with a 5-mm-long Silastic capsule (inner diameter, 
1.57 mm; outer diameter, 2.41 mm) filled with testoster- 
one at around the time of puberty (PN45). This design 
allows us to attribute any differences between XX and 
XX + T mice only to the neonatal hormone exposure. 
Post-weaning, animals belonging to the same litter were 
housed in the same cage (maximum of three animals per 
cage), although males and females were separated. As all 
pups in a litter received the same treatment (either T or 
vehicle), all adult animals housed in the same cage had 
undergone the same hormonal manipulations. 

We had a total of 12 experimental groups: XX, XY, 
and XX + T tissue from two ages and two distinct brain 
regions. Each experimental group had two biological 
replicates, each comprising a pool of tissue from three 
animals, making for a total of 24 samples. 

Brain dissections 

At PN4 or PN60, whole brain was rapidly removed from 
the skull and brain regions of interest were dissected 
under a microscope on an ice-cold slide. After removal 
of the dura mater, two cuts through the brain along the 
coronal plane were made. The first was at the midpoint 
of the optic chiasm (0.14 mm anterior to bregma) and 
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the second was where the optic tract enters the brain 
(0.58 mm posterior to bregma). The resulting slab of tis- 
sue was then placed posterior side down. The BNST/ 
POA was defined as the region ventral to the lateral ven- 
tricle and bounded laterally by the medial edge of the in- 
ternal capsule. The striatum was defined as the tissue 
between the external capsule and the anterior commis- 
sure, bounded laterally by the cortex and medially by the 
internal capsule. After dissection, the tissue was immedi- 
ately placed on dry ice and stored at -80°C until it was 
processed for downstream experiments. 

Measurement of levator ani and bulbocavernosus muscles 

The dose of testosterone used in this study has previ- 
ously been shown to completely masculinize BNST vol- 
ume and cell number [24]. To confirm the efficacy of 
hormone treatments here, we randomly selected a subset 
of the XX (N= 4) and XX + T (N= 5) mice on PN4 and 
examined the size of the androgen-sensitive levator ani 
and bulbocavernosus muscles. The perineums were 
processed as described previously [25], and the maximal 
cross-sectional area of the levator ani was determined 
by tracing around the muscle using Stereolnvestigator 
(MBF Biosciences, Williston, VT, USA) software. 

Reduced representation bisulfite sequencing library 
construction 

Genomic DNA from mouse brains was extracted for 
making reduced representation bisulfite sequencing 
(RRBS) libraries following the standard RRBS protocol 
[26]. The genome was digested with the Mspl enzyme, 
a methylation-insensitive restriction enzyme. These 
Mspl-digested samples were ligated with Illumina 
adaptors (Illumina Inc., San Diego, CA, USA) and size- 
selected. Fragments from 100 to 200 bases were se- 
lected as these are enriched for CpG-rich regions, such 
as CpG islands, promoter regions, and enhancer ele- 
ments. In total we selected 500 K distinct fragments 
for sequencing. These fragments were denatured and 
treated with sodium bisulfite to reveal their methyla- 
tion status (CpGenome Universal DNA Modification 
Kit, Cat. No. S7820, Millipore, Billerica, MA, USA). 
Libraries were then polymerase chain reaction (PCR)-amp- 
lified with MyTaqHS (BIO-25047, Bioline, Taunton, MA, 
USA) and sequenced using the Solexa sequencing technol- 
ogy (Illumina Hiseq 2000 sequencers). The reads were 
aligned to the reference genome (mouse mm9) using the 
modified bisulfite aligner, BS Seeker [27]. To generate 
genome-wide DNA methylation profiles, we calculated the 
methylation level for each covered cytosine on the genome. 
As bisulfite treatment converted unmethylated cytosines 
(Cs) to thymines (Ts), we estimated the methylation level at 
each cytosine by #C/(#C + #T), where #C is the number of 
methylated reads and #T is the number of unmethylated 



reads. This number represents the average methylation 
level at that particular site across the cell population tested. 
In this study we only included cytosines that were covered 
by at least four reads for the analysis. 

Identifying differentially methylated regions and 
differentially methylated genes 

We first searched for differentially methylated regions 
(DMRs) that showed significant differential methylation. 
We grouped nearby cytosine sites into units called frag- 
ments. In each pairwise comparison, a Student's t test 
was performed at each site to quantify the difference be- 
tween the groups at that site. This generated a £-score 
that represented the difference between the groups (the 
larger the ^-scores, the more different the methylation 
levels in that pairwise comparison). In order to get an 
accurate measurement of these differences after the sites 
were combined into fragments, the ^-scores of all sites in 
that fragment are averaged to produce a z-score. To 
qualify as DMRs, the fragment had to (1) show a differ- 
ence of >10% in mean methylation level between the 
two groups being compared, (2) have at least three cyto- 
sines for which methylation levels were observed in all 
relevant samples, and (3) have a z-score below a thresh- 
old relevant to that comparison. The selection of the 
z-score threshold was based on the false discovery rate 
estimated by comparing the real data to simulated 
methylomes as the control for false discovery rate (FDR) 
computation (full procedure below). These DMRs were 
then associated with a gene if there was a transcription 
start site within 5 kb of them or if they overlapped with any 
known genes to identify differentially methylated genes. We 
used GeneVenn (http://genevenn.sourceforge.net/, [28]) to 
determine overlap between gene sets. 

Estimating false discovery rate 

To assess the false discovery rate for our DMRs, we con- 
structed simulated methylomes, with the same read 
coverage per site as the real samples. For each CG site 
in each simulated sample, we then simulated the reads 
(C if methylated or T if unmethylated) based on the 
average methylation level (Pm) from all real samples at 
this CG site. The number of methylated reads (Cs) at a 
site of coverage n is a random sample from the binomial 
distribution B(n, Pm). We repeated our simulation of 
reads throughout the genome for all samples. The result- 
ing samples had the same average methylation levels as 
the real sample. Since the reads were simulated from the 
binomial distribution with the same average methylation 
levels as in the real samples, the differences in methyla- 
tion patterns across genes, repeats, promoters, etc. were 
preserved. The simulated data has the same coverage as 
the real samples so the statistical power is not affected. 
The simulated methylomes should have no difference in 
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methylation levels between the two comparison groups 
(i.e., no DMRs), since they are all selected using the 
same methylation frequency. Any DMRs (and the DMR- 
associated genes) identified from these simulated sam- 
ples are thus considered false positives. Then, for each 
comparison we repeated the whole procedure to detect 
the DMR on simulated samples: we first performed t tests 
on individual sites and then summarized the ^-scores per 
fragment with a z-score. For each z-score threshold, we 
computed the numbers of DMRs that were found in the 
simulated data to those found in the real data. We used the 
ratio of these to compute the FDR. We chose a z-score 
threshold that resulted in a false discovery rate less than 
10% in all comparisons. 

Traditional (Sanger) bisulfite sequencing 

Strand- specific primers were designed for the bisulfite- 
converted genome of the region of interest and synthe- 
sized by Life Technologies (Carlsbad, CA, USA). For 
Micalll (target region, chrl5:78965961-78966079; nega- 
tive strand, mm9), the forward primer was ATTTTTGT 
TATTGGGAAGGATAAGG and the reverse primer was 
AAACCCCAACCATACATAATCTCTA. The cycling con- 
ditions were 1 cycle at 95°C for 2 min; 35 cycles at 95°C for 
30 s, 58°C for 1 min, and 60°C for 1 min; and 1 cycle at 
60°C for 15 min. For Fzd9 (target region, chr5: 135725414- 
135725524; negative strand, mm9), the forward primer was 
TGAATTGATTGGGTTTTGTTATGTA and the reverse 
primer was ACTAATAATACCCACCACCAAAAAC. The 
cycling conditions were 95°C for 2 min; 40 cycles at 95°C 
for 30 s, and 60°C for 2 min; 1 cycle at 60°C for 15 min. 
The samples used for traditional bisulfite sequencing (n = 
2-3 per experimental group) were not used in the RRBS 
experiments although they were generated at the same 
time. One microgram for each sample was bisulfite-treated 
and purified (CpGenome Universal DNA Modification Kit, 
Cat. No. S7820, Millipore). Forty nanograms of bisulfite- 
converted DNA was used in each PCR (MyTaqHS, BIO- 
25047, Bioline). After gel purification, amplicons were 
cloned into pCR4-TOPO TA vectors (TOPO TA Cloning 
Kit for Sequencing, K4575-01SC, Life Technologies). 
Fifteen to twenty colonies were sent for Sanger sequen- 
cing using the M13R primer (Laragen Inc., Culver City, 
CA, USA). 

Heat maps of methylation levels in differentially 
methylated regions 

A union set of DMRs was collected from all pairwise 
comparisons (all ages, sex, and regions). Among them, 
we selected 4,086 DMRs that were extremely differentially 
methylated (longer than 50 bp, delta methylation >25%, 
and -z-score >3.5). These are potentially regions that 
are susceptible to methylation changes. The average 
methylation levels in all XX, XX + T, and XY groups 



were plotted in heat maps with hierarchical clustering 
of the DMRs. 

Gene ontology using Ingenuity Pathway Analysis 

Functional analysis of statistically significant CpG 
methylation changes was performed with Ingenuity 
Pathway Analysis (IPA; Ingenuity Systems, www.ingenu- 
ity.com). Ingenuity functional analysis identified net- 
works, canonical signaling pathways, and biological 
functions and/or diseases that were most significantly af- 
fected by testosterone and age. For all analyses, datasets 
containing gene identifiers and corresponding delta 
methylation values were uploaded into IPA. The genes 
were overlaid onto a global molecular network devel- 
oped from information in the Ingenuity Pathway Know- 
ledge Base. Networks of these focus genes were then 
algorithmically generated based on their connectivity. 
To identify biological functions and diseases that were 
enriched in the different datasets, genes were associated 
with biological functions and/or diseases in the Ingenu- 
ity Knowledge Base. Right-tailed Fisher's exact test was 
used to calculate a p value determining the probability 
that each biological function and/or disease assigned to 
that dataset was due to chance alone. In this method, 
the p value for a given process is calculated by investi- 
gating (1) the number of participating genes in that 
process and (2) the total number of genes known to be 
related to that process in the selected reference set. The 
more genes involved, the more significant the p value. 
Canonical pathways analysis identified the pathways 
from the IPA library of pathways that were most signifi- 
cant to the dataset. The significance of the association 
between the dataset and the canonical pathway was 
measured in two ways: (1) a ratio of the number of mol- 
ecules from the dataset that map to the pathway divided 
by the total number of molecules that map to the canon- 
ical pathway was determined and displayed in the tables 
that follow; (2) Fisher's exact test was used to calculate a 
p value determining the probability that the association 
between the genes in the dataset and the canonical path- 
way was explained by chance alone. 

Testosterone measurements 

Samples were collected at the time of euthanasia. In all 
cases, blood was obtained from the carotid artery follow- 
ing decapitation. Blood samples were then processed 
to isolate serum and stored at -20°C until assays for 
testosterone were performed. Testosterone assays using 
radioimmunoassay were performed by Ligand Assay and 
Analysis Core at the University of Virginia Center for Re- 
search in Reproduction (supported by NICHD (SCCPIR) 
Grant U54-HD28934). Testosterone measurements were 
performed in singlet reactions using Siemens Medical Solu- 
tions Diagnostics testosterone RIA (Siemens Healthcare, 
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Malvern, PA, USA) with a reportable range of 47.3- 
170.5 ng/L. There were no significant differences in the 
measured testosterone levels between our experimental 
groups using the Kruskal-Wallis one-way analysis of vari- 
ance test (H= 3.8, 2df,p = 0.15). 

Processing for gene expression analysis 

Total RNA samples were derived from five pools of three 
animals using AllPrep DNA/RNA Micro Kit (Qiagen, 
Valencia, CA, USA) according to the manufacturer's 
protocol. This kit enabled simultaneous purification of 
genomic DNA and total RNA from all tissues. Samples 
were prepared from the same tissue samples that were 
used to create the RRBS libraries. RNA was quantified 
using a Ribogreen fluorescent assay (Life Technologies) 
and normalized to 10 ng/|il prior to amplification. 
Amplified and labeled cRNA was produced using the 
Illumina specific Ambion TotalPrep kit (Applied Bio- 
systems Inc., Foster City, CA, USA). First and second 
strand cDNA were produced using the Ambion kit and 
purified using a robotic-assisted magnetic capture step. 
Biotinylated cRNA was produced from the cDNA tem- 
plate in a reverse transcription reaction. Typical yields 
were in excess of 1.5 |ig. After a second Ribogreen 
quant and normalization step, amplified and labeled 
cRNA was hybridized overnight at 58°C to MouseRef-8 
vl.l BeadChip expression arrays from Illumina. To 
minimize array-to-array variability, a cRNA sample 
from each of the experimental groups was hybridized 
to each of the beadchips (n = 5/group) according to the 
manufacturer's protocol. The MouseRef-8 vl.l beadchip 
contains over 24,000 well-annotated RefSeq transcripts and 
allows eight samples to be interrogated in parallel. 
Hybridization was followed by washing, blocking, staining, 
and drying on the Little Dipper processor (Agilent 
Technologies, Sta. Clara, CA, USA). Array chips were 
scanned on either the BeadArray reader or the iScan 
reader (Illumina). 

Microarray data analysis of gene expression 

cDNA microarray data analysis was performed using the R 
software and Bioconductor packages. The raw intensity 
data were first log 2 transformed and then the outlier 
samples were detected based on the average inter-sample 
correlations. Three samples with average inter-sample cor- 
relations larger than two standard deviations (SDs) from 
the mean across all samples were removed from follow-up 
analysis. Then the samples were normalized by quantile 
normalization. Probes were considered robustly expressed 
if the detection p value was <0.05 in at least 20 samples 
(the total number of samples for every genetic/treatment 
group) in the dataset. After the data quality controls, 57 
samples with 13,776 probes were retained for differential 
expression analysis. Hierarchical clustering analysis using 



one minus inter-array correlation measures was further 
performed to assess the sample clustering patterns. 

Differential expression analysis was conducted using 
the R limma package [29], and unless otherwise speci- 
fied, the significance threshold was Benjamini-Hochberg 
(BH) -adjusted p value <0.05. Limma uses linear models 
to robustly assess the differential expression of genes. 
Histograms of p values were plotted to show if there was 
a significant differential expression signal genome-wide. 
To identify the genes that differed by age in the different 
groups/brain regions, the factorial design was applied 
using limma linear models. The genes were identified as 
showing significant dynamic changes with age if their 
p values associated with the interaction term between 
the factors of group and age were less than 0.005. 

Quantitative reverse transcription-polymerase chain reaction 

Reverse transcription was performed using the Tetro 
cDNA Synthesis Kit (Bioline, catalog no. BIO-65043) 
with 1 \ig of total RNA as template. The RNA samples 
used for validation were from the original microarray 
samples contingent on availability (n = 3-4 per geno- 
type). The sequences of PCR primers are as follows: for 
Alcam, forward primer 5'-CGA ACC CTG CCT GTG 
TCA TGC ACA ATA-3 ' and reverse primer 5'-TAT 
CGT CTG CCT CAT CGT GCT CTG GAA T-3'; for 
Gapdh, forward primer 5-TGC CGC CTG GAG AAA 
CC-3' and reverse primer 5 -CCC TCA GAT GCC 
TGC TTC AC-3'; for Gatadl, forward primer 5'-GAA 
ATT CAC AG A AGG TCG GC-3' and reverse primer 

5 -AAT ATA CTC CCT TGT AGA AGA TTG A-3'; for 
Utx, forward primer 5' -CCA ATC CCC GCA GAG 
CTT ACC T-3' and reverse primer 5 -TTG CTC GGA 
GCT GTT CCA AGT G-3'. All primers used spanned at 
least one intron. Glyceraldehyde-3-phosphate dehydrogen- 
ase (GAPDH) was used for normalization of gene expres- 
sion between samples. Quantitative reverse transcription- 
polymerase chain reactions (qRT-PCRs) were carried out in 
duplicate utilizing the Sybr Green-based SensiFAST™ SYBR 

6 Fluorescein Kit (Bioline, catalog no. BIO-96005). The 
thermocycling parameters for all reactions were 1 cycle at 
95°C for 2 min, 40 cycles at 95°C for 10 s, 60°C for 10 s, 
and 72°C for 10 s. We used the standard curve method to 
determine relative expression and assessed significance 
using the Student's t test (a = 0.05). Data are expressed as 
fold change where the expression level in the XX group has 
been set to 1. 

Data access 

The data discussed in this publication have been deposited 
in NCBI's Gene Expression Omnibus and are accessible 
through GEO Series accession number GSE50218 (http:// 
www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE50218). 
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The expression data generated for this study have been 
submitted to the NCBI Gene Expression Omnibus (GEO; 
http://www.ncbi.nlm.nih.gov/geo/) under accession number 
GSE49986. 

Results 

Sex effects on genome-wide methylation data 

To examine sex-specific and hormone-induced changes in 
brain DNA methylation, we compared genome-wide maps 
of DNA methylation in adult mouse striatum and BNST/ 
POA between XX and XY mice generated using RRBS. All 
animals were gonadectomized and provided with identical 
hormone replacement prior to sacrifice so that any differ- 
ences seen could not be due to differences in circulating 
hormones at the time of sacrifice (see 'Methods'). We had a 
total of 12 experimental groups: XX, XY, and XX + T tissue 
from two ages and two distinct brain regions. Each experi- 
mental group had two biological replicates, each compris- 
ing a pool of tissue from three animals, making for a total 
of 24 samples. Forty-six percent of our sequencing reads 



uniquely mapped back to the mouse genome, which re- 
sulted in an average sequencing depth of 58 x at each CpG. 
We assayed approximately 1.39 million CpGs in 125 mil- 
lion uniquely mapped reads in each sample (Additional file 
1: Table SI). The majority of samples had >1 million dis- 
tinct CpGs covered with a minimum of 4x coverage. CpG 
sites that were not present in all comparison groups were 
excluded from further analysis. 

We first examined global levels of CpG, CHH, and 
CHG methylation. At CHH sites, any two bases can fol- 
low the initial C, whereas at CHG sites, any base can be 
in the second position but the third base must be a G. 
There were no significant effects on these global levels 
due to sex or hormonal treatment. However, the age of 
the animal was positively correlated with methylation 
level, especially at non-CG sites (CHG p = 9.19E - 8 in 
the striatum, p = 1.37E - 8 in the BNST/POA; CHH p = 
4.88E - 10 in the striatum, p = 6.88E - 11 in the BNST/ 
POA; Student's t test for all comparisons), which is con- 
sistent with recent studies [12,13] (Figure la,b,c). Global 
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Figure 1 DNA methylation patterns in XX, XY, and XX + T. Non-CpG methylation increases with age. (a) Global CpG methylation levels are 
significantly higher at PN60 than PN4 in the BNST/POA but not the striatum. CHG (b) and CHH (c) methylation levels are significantly higher at 
PN60 than at PN4 in both brain regions. (d, e) The heat maps show correlations between the methylation levels across all samples in the (d) 
striatum and (e) BNST/POA. The diagonal axis running from the bottom left to top right comer is the line of symmetry where each sample correlates with 
itself. The brightest yellow shade represents full correlation (=1) whereas the darkest red shade represents the lowest correlation value (0.98; key at top left). 
Complete clustering is based on the Euclidean distance, (f-i) These heat maps show differentially methylated regions (DMRs) driven by sex and/or 
testosterone at PN60. The color key shows absolute methylation level from 0% (blue) to 100% {red). The azure trace represents the moving average 
of the methylation levels in the neighboring windows, (f) Autosomal DMRs in the striatum, (g) Autosomal DMRs in the BNST/POA. (h) X-linked 
DMRs in the striatum, (i) X-linked DMRs in the BNST/POA. 
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CpG methylation levels were significantly increased with 
age in the BNST/POA (p = 5.9E - 3 by the Student's t test) 
but not the striatum. 

We then focused our analysis on CpG methylation. 
Overall, the CpG methylation profiles of adult XX and 
XY striatum and BNST/POA were highly similar across 
all chromosomes (Pearson coefficient, 0.99) (Figure ld,e), 
indicating that the genomic profiles of 5-mC are highly 
reproducible. Despite this overall similarity, hierarchical 
clustering clearly identified sex- and testosterone-driven 
(Figure lf,g,h,i) methylation differences. Testosterone- 
driven differences (different in XX vs. XX + T) were 
more abundant on autosomes (Figure lf,g). In many re- 
gions showing these differences, XX + T CpG methyla- 
tion patterns resemble those of XY. On the other hand, 
sex-specific differences (XX vs. XY) were clearest on the 
X chromosome, as expected, given the involvement of 
CpG methylation in the process of X inactivation and 
XX + T tended to resemble XX (Figure lh,i). 

To identify genes that undergo sex-specific methyla- 
tion, we compared the methylomes of male (XY) and fe- 
male (XX) mice at PN60 and found a large number of 
genes that showed sex differences in methylation pat- 
terns in both the striatum (1,577 genes) and the BNST/ 
POA (1,026 genes) (Additional file 2: Table S2). Of the 
1,577 genes in the striatum that showed sex-specific methy- 
lation, 420 were on the X chromosome and 348 (82.86%) of 
these were more methylated in females (Figure 2a). 
Similarly, in the BNST/POA, 426 of the 1,026 sex- 
affected genes were X-linked and 359 (84.27%) of these 
demonstrated elevated methylation levels in females 
(Figure 2b). This is likely due to X chromosome inacti- 
vation in females [30,31]. Focusing on the effect of sex 
on autosomes, we identified 1,157 and 600 genes with 
sex differences in the adult striatum and BNST/POA, 
respectively. Interestingly, a substantial number of these 
genes showed higher methylation in XY than XX. In the 
striatum this proportion was 1,070/1,157 (92.48%) of 
autosomal genes (Figure 2a). In the BNST/POA, 520/ 
600 (86.67%) genes were more methylated in XY relative 
to XX during adulthood (Figure 2b). We next compared 
the set of genes that showed sex-specific methylation in 
each brain region to each other. Only genes that showed 
sex-specific methylation differences at the same chromo- 
somal location and in the same direction were consid- 
ered to be common to both brain regions. Based on 
these criteria, we found 362 genes, which we term dual 
region sex-specific genes (Figure 2c). There was a sub- 
stantial overrepresentation of X-linked genes (286/362 
or 79%). Together, these data show that methylation sta- 
tus of many genes in the brain is subject to sex- and 
region-specific regulation and that the process of X 
inactivation largely affects the same genes in both the 
striatum and the BNST/POA. 



a 



Striatum 



BNST/POA 





■ X t Me in XX 
X t Me in XY 
Autosome f Me in XX 

■ Autosome f Me in XY 
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Dual region sex-specific genes 
PN60 



1215 



362 664 



Striatum^ 


BNST/POA 


(XX vs. XY) 


(XX vs. XY) 


21% \_ 79% 



Autosomal X-linked 

Figure 2 CpG methylation at many genes is sexually dimorphic. 

Displayed are the fractions of X and autosomal genes displaying 
sexually dimorphic CpG methylation in (a) PN60 striatum and (b) 
PN60 BNST/POA. (c) Genes with sexually dimorphic CpG methylation 
at PN60 that are common to both the striatum and BNST/POA. 
The Venn diagram shows the overlap {green) between sexually 
dimorphic genes from the striatum {blue) and BNST/POA {yellow). The 
numbers of genes belonging to each category are shown within the 
relevant circles. The bar below the Venn diagram shows the percentage 
of autosomal vs. X-linked genes in the overlap region. 



Testosterone-induced modification of brain CpG methylation 

Sex differences in methylation patterns can be attributed 
to differences in sex hormones or to other influences 
(e.g., direct genetic effects of the sex chromosomes). To 
determine the effects of neonatal testosterone exposure 
on CpG methylation genome-wide, we gave newborn fe- 
male mice a masculinizing dose of testosterone (XX + T) 
and compared their CpG methylation patterns to that of 
control females (XX). The masculinizing effect of this 
dose of testosterone was verified by examining the le- 
vator ani and bulbocavernosus muscles. These muscles 
are androgen-sensitive and while they are present in 
both sexes at birth [25], they die postnatally in females 
due to the lack of testosterone. The maximal cross - 
sectional area of the levator ani was measured in PN4 
animals. Testosterone treatment resulted in a clear spar- 
ing of this muscle with it being five times larger in XX + 
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T than in XX (269,150 ± 11,908 jim 2 in XX + T vs. 
53,721 ± 9,414 (W in XX; t= 13.6, p <0.0001) (Figure 3a; 
Additional file 3: Table S3). Additionally, the bulbocaver- 
nosus muscle was absent in all XX animals examined 
but preserved in all XX + T. 

We found that neonatal testosterone exposure altered 
the methylation status of a substantial number of gen- 
omic fragments. These fragments mapped to a relatively 
small number of genes at PN4 (68 genes in the striatum 
and 45 genes in the BNST/POA) (Figure 3b). However, 
by day 60 of life, a much larger number of genes showed 
methylation changes in response to neonatal testoster- 
one (1,377 and 740 genes differed in methylation status 
in testosterone-treated PN60 females relative to control 
PN60 females in striatum and BNST/POA, respectively) 
(Figure 3b; Additional file 3: Table S3). A substantial 
fraction of these genes displayed increased 5-mC in re- 
sponse to testosterone. In the striatum, 51/68 (75.00%) 
genes at PN4 and 1,324/1,377 (96.15%) genes at PN60 
showed greater methylation in XX + T relative to XX 
mice (Figure 3c). Similarly, in the BNST/POA 38/45 
(84.44%) genes at PN4 and 705/740 (95.27%) genes at 
PN60 also showed greater methylation in female mice 
treated with testosterone (Figure 3c; Additional file 3: 
Table S3). 

Of the testosterone-affected genes, 654/1,377 (47.49%) 
and 265/740 (35.81%) in the striatum and BNST/POA, 
respectively, also exhibited male vs. female differences 
in their methylation patterns at PN60 (Figure 3d,e; 
p value = 1.00E - 111, hypergeometric test; fold enrich- 
ment over chance, 4.55 (striatum), 6.07 (BNST/POA); 
Additional file 4: Table S4). There are several possible 
reasons why many of the genes whose methylation 
status was changed by testosterone did not show sex 
differences: (1) a lack of statistical power such that 
some differences are not confirmed in two different 
comparisons, (2) a pharmacological effect of the rela- 
tively high dose of testosterone used here, or (3) testos- 
terone effects found in XX but not normally in XY 
mice. For some dependent measures, the presence of a 
Y chromosome has effects opposite to that of testoster- 
one [32]. In the current context, this means that some 
effects of early exposure to androgens on CpG methy- 
lation in males may be masked by Y chromosome 
factors. 

To validate our RRBS findings, we performed trad- 
itional (Sanger) bisulfite sequencing. We chose one 
testosterone-affected locus from each brain region for 
validation. For the striatum, we chose Micalll (higher 
methylation in XX than XX + T). For the BNST/POA, 
we chose Fzd9 (higher methylation in XX + T than 
XX). The results from traditional bisulfite sequencing 
were in agreement with the RRBS data (Additional 
file 5: Figure SI). 



Brain region specificity of testosterone's effects 

Next, we compared the testosterone-affected genes in 
the striatum to those in the BNST/POA. As with our 
identification of dual region sex-specific genes, only 
genes where testosterone-affected methylation at the 
same chromosomal location and in the same direction 
were considered. At PN4, only two genes (Sorcs2 and 
LonrfS) were affected by testosterone in both the stri- 
atum and the BNST/POA (Figure 3f). At PN60, this 
number rose to 92, (neither Sorc2 nor LonrfB were 
present in this list) (Figure 3g). Unlike the list of PN60 
dual region sex-specific genes, the X chromosome was 
not substantially over-represented in this comparison 
(16 of the 92 dual region testosterone-affected genes 
were X-linked). This result may be expected since in 
both XX and XX + T animals, one X chromosome is ex- 
pected to undergo inactivation and the same genes are 
inactivated across tissues [33]. These data indicate that 
testosterone's effects are more region-specific than that 
of sex (as there are only 92 dual region testosterone- 
affected genes compared to 362 sex-specific ones) and 
that genetic sex has a much stronger effect on methyla- 
tion of X-linked genes, whether they were subject to or 
escaped X inactivation (compare Figure 2c to 3g). 

In both the striatum and BNST/POA, the sexually 
dimorphic gene sets showed a greater proportion of X- 
linked genes compared to the testosterone-affected gene 
sets (PN60 striatum, 26.6% vs. 6.1%; PN60 BNST/POA, 
41.5% vs. 9.7%). When we consider just inactivation es- 
capees, the results are largely similar to those from the 
whole set. At both ages and in both brain regions, we 
detected 5 of 14 known X inactivation escapees as being 
sexually dimorphic in their methylation (1810030O07Rik, 
6720401G13Rik, Bgn, Midi, Shroom4) [34,35]. In contrast, 
only 1 of these 14 (striatum, Midi; BNST/POA, 
672040 lG13Rik) were affected by testosterone and 
only at PN60. Thus, it appears that the effect of testos- 
terone on X-linked genes is rather limited in both the 
striatum and BNST/POA. 

To examine the characteristics of testosterone-affected 
genes in adulthood, we used Ingenuity Pathway Analysis 
(IPA, Ingenuity Systems). Our analysis revealed that tes- 
tosterone alters the methylation of genes belonging to a 
wide range of biological processes and functions (Table 1; 
Additional file 6: Table S5). Functional categories related 
to nervous system development were strongly repre- 
sented in the datasets from both the striatum and the 
BNST/POA. Many of these functional categories were 
common to both regions of the brain and represented 
basic processes that are crucial for general neural func- 
tion (for example, morphology of nervous tissue, neuro- 
genesis, guidance of axons, and morphology of dendritic 
spines). Genes related to organization of cytoskeleton, 
microtubule dynamics, and apoptosis were also altered 



Ghahramani et al. Biology of Sex Differences 2014, 5:8 
http://www.bsd-journal.eom/content/5/1/8 



Page 9 of 1 8 



a 




XX + T 



1600 
1400 
1200 
1000 
800 
600 
400 
200 
0 



1377 






740 




J 






PN4 PN60 


PN4 PN60 


Striatum 


BNST/POA 


Striatum 







e 



BNST/POA 




XX vs. XY 



XX vs. XX + T XX vs. XY 



XX vs. XX + T 



f g 

Dual region testosterone-affected genes Dual region testosterone-affected genes 



PN4 




PN60 




Striatum 
(XX vs. XX + T) 



Figure 3 Perinatal testosterone treatment affects CpG methylation at numerous loci, (a) The dose of testosterone given to female mice at 
PNO is masculinizing. Representative pictures of the levator ani (LA) are shown for PN4 XX {left panel) and XX + T {right panel). The muscle size is 
reduced in the absence of testosterone. The arrows indicate the outline of the LA, which has a much larger surface area in XX + T than in XX 
(269,1 50 ± 1 1,908 urn 2 in XX + T vs. 53,721 ±9,414 urn 2 in XX; f= 13.6, p <0.0001). (b) Number of genes where methylation is altered in the XX 
vs. XX + T comparison in PN4 and PN60 striatum and PN4 and PN60 BNST/POA. (c) Percentage of genes that exhibit testosterone-dependent hypo- or 
hyper-methylation at each age in each tissue. (d, e) The Venn diagrams show the overlap {green) between genes with sexually dimorphic methylation 
patterns {blue) and those that are affected by testosterone {yellow) in the (d) striatum and (e) BNST/POA. The numbers of genes belonging to each 
category are shown within the relevant circles. (f, g) The number of genes that are affected by testosterone in both the striatum and the BNST/POA at (f) 
PN4 and (g) PN60. The overlap between testosterone-affected genes from the striatum {blue) and BNST/POA {yellow) is shown in green. The numbers of 
genes belonging to each category are shown within the relevant circles. The bar below the Venn diagram shows the percentage of autosomal vs. X-linked 
genes in the overlap region. 
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Table 1 Biological functions associated with differentially methylated genes in XX vs. XX + T mice determined by IPA 


Category 


Function annotation 




Striatum 




BNST/POA 






p value 


Number of genes 


p value 


Number of genes 


Nervous system development 
and function 


Morphology of nervous system 
Development of central nervous system 


1 .02E - 08 
6.63E-08 


108 
80 


3.91 E- 06 
4.71 E- 08 


63 
54 




Morphology of nervous tissue 


7.88E - 06 


72 


7.75E-05 


44 




Neuritogenesis 


1.05E-05 


48 


1.03E-07 


37 




Outgrowth of neurites 


4.68E - 05 


39 


2.19E-03 


22 




Coordination 


4.69E - 05 


31 


1 .67E - 02 


15 




Axonogenesis 


1 .48E - 04 


24 


1.18E- 04 


17 




Excitatory postsynaptic potential 


1 .5 1 E - 04 


18 


5.55E-03 


10 




Growth of neurites 


1.75E-04 


41 


5.21 E- 03 


23 




Morphology of neurites 


7.77E - 04 


23 


1 .99E - 04 


17 




Morphology of dendritic spines 


7.78E - 04 


5 


7.85E-04 


4 




Guidance of axons 


9.29E - 04 


21 


1 .60E - 05 


18 




Outgrowth of axons 


9.56E - 04 


14 


9.82E - 03 


8 




Neurotransmission 


6.32E-05 


44 


_ 


_ 




NMDA-mediated synaptic current 


2.93E-04 


6 


- 


- 




Action potential of cells 


o cop n/i 


1 Q 








Long-term depression 


1 30E — 03 


1 3 








Quantity of neurons 






3.81 E — 05 


24 




Cell viability of neurons 






1.93E-02 


12 


L.cllUldl dbbclllUiy dl IU UlLjdl IIZdLIUI 1 


wiydl IIZdLIUI 1 Ul Ly LUbKclcLUI 1 


I .ooc u/ 


1 9^ 
I ZD 




77 
/ / 




Microtubule dynamics 


2.17E-07 


108 


7.64E - 06 


65 


Skeletal and muscular system 
development and function 


Development of muscle 


2.55E-06 


49 






Cell death and survival 


Apoptosis 


7.68E - 06 


250 


1 .70E - 03 


138 


Cellular growth and proliferation 


Proliferation of neuronal cells 






1.25E-02 


13 


Behavior 


Learning 
Social behavior 


3.33E-03 


38 


2.81 E- 04 
2.59E-03 


27 
7 



This table shows examples of functional 
XX + T comparison for both the striatum 



categories enriched in the XX vs. 
and BNST/POA. 



in adulthood by neonatal testosterone treatment in both 
brain regions. We observed some biologically relevant 
differences between the striatum and the BNST/POA. In 
the XX vs. XX + T comparison, functional categories 
uniquely enriched in the adult striatal dataset included 
neurotransmission, NMDA-mediated synaptic current, 
action potential of cells, long-term depression, and devel- 
opment of muscle. The striatum also showed an enrich- 
ment of functional categories related to neurological 
disease (e.g., dyskinesia, Huntington's disease, hyperactive 
behavior, seizures, mood disorder, ataxia, and amyo- 
trophic lateral sclerosis) (Table 2). Functional categories 
unique to the BNST/POA included quantity of neurons, 
cell viability of neurons, and proliferation of neuronal 
cells (Table 1). The presence of genes related to expan- 
sion and survival of neurons in the BNST/POA 
testosterone-affected dataset is intriguing given the role 



that cell death plays in the establishment of sex differ- 
ences in this region [22] and recent evidence that cell 
birth may be an important mechanism that helps main- 
tain the sexual dimorphism in the sexually dimorphic 
nucleus of the POA [36]. One effect of this cell death is 
that cells sampled at PN4 in XX are most likely not the 
same population as the cells sampled in adulthood. 
Therefore, the methylation differences that we observe 
might be a result of this change in the cellular popula- 
tion between the two time points studied. 

We next examined the genomic features of fragments 
identified by RRBS as differing between XX and XX + T 
at PN60. In both brain regions, we found an underrepre- 
sentation of methylation changes not only in promoters 
(defined as ±500 bp relative to the transcription start 
site) but also in CpG islands. Gene-body methylation 
(the entire gene from the transcription start site to the 
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Table 2 Examples of top 'neurological disease' functional 
categories that were significantly enriched in striatal 
testosterone-influenced genes 

Function annotation p value Number of genes 

Movement disorder 1.03E-04 118 

Congenital anomaly of brain 6.81 E- 04 26 

Dyskinesia 3.36E-03 73 

Huntington's disease 3.43E-03 70 

Hyperactive behavior 3.47E-03 18 

Schizophrenia 3.90E-03 56 

Jervell and Lange-Nielsen syndrome 4.66E-03 2 

Seizures 5.08E-03 33 

Mood disorder 5.35E-03 50 

Hydrocephalus 6.76E-03 11 

Ataxia 6.89E - 03 22 

Amyotrophic lateral sclerosis 7.00E-03 18 

Incoordination 9.03E-03 3 

Oligodendroglioma 9.64E-03 5 

Spina bifida 1.33E-02 7 

Degeneration of brain 1.34E-02 8 

Examples of top 'neurological disease' functional categories that were 
significantly enriched in the list of testosterone-influenced genes (XX vs. 
XX + T) in the striatum as determined by Ingenuity Pathway Analysis. 

end of the transcript) contributed substantially to testoster- 
one-altered CpGs. CpGs located within introns were most 
susceptible to changes by testosterone (p value <L00E - 94; 
Fisher's exact test). Testosterone-modified CpGs were also 
over-represented in exonic regions (p value <L00E -17; 
Fisher's exact test). These results are consistent with recent 
reports which suggest that in mammals, cell type-specific 
CpG methylation-related gene regulation mostly occurs at 
alternative promoters within gene bodies [37]. When we 
compared the methylation patterns of three subclasses 
of repetitive elements— short interspersed elements 
(SINEs), long interspersed elements (LINEs), and simple 
repeat regions— in our testosterone dataset, we detec- 
ted a depletion of testosterone-affected CpG sites in 
LINEs {p value <1.00E- 14; Fisher's exact test). 

Identification of stably differentially methylated genes 

To identify genes that showed stable differential methy- 
lation due to testosterone, we compared the genes that 
showed differential methylation between XX vs. XX + T 
at PN4 to the XX vs. XX + T genes at PN60. We defined 
stably differentially methylated genes as those in which a 
difference in CpG methylation at PN4 was also found at 
PN60. We found a limited number of such genes (19 
genes in the striatum and 11 genes in the BNST/POA) 
(Table 3). About a third of these genes are X-linked but 
none are known to escape X inactivation, implying that 
testosterone does not affect X inactivation status. In 14/ 



19 genes in the striatum and 10/11 genes in the BNST/ 
POA, the magnitude of testosterone-driven changes in 
methylation (delta methylation between XX and XX + T) 
was similar at both PN4 and PN60. Except for a few 
genes (three genes in the striatum and two genes in 
the BNST/POA), the direction of the testosterone- 
dependent epigenetic modifications remained the same 
between PN4 and PN60. The majority of the locations 
where stable testosterone-driven differential methylation 
is observed are associated with multiple genomic features. 
In both brain regions, the genomic feature most strongly 
associated with stable methylation changes is gene bodies 
(associated with 18/19 genes in the striatum and 11/11 
genes in the BNST/POA). Next is methylation changes in 
CpG islands (striatum, 11/19 genes, BNST/POA, 6/11 
genes) and then promoter regions (striatum, 7/19 genes, 
BNST/POA, 6/11 genes). Overall, the small number of 
genes showing stable differential methylation indicates 
that testosterone-induced methylation patterns are dyna- 
mic and can develop long after the initial exposure. 

Testosterone-induced masculinization of methylation 

We next assessed whether testosterone induces a broad 
shift in CpG methylation in the brains of XX + T mice 
from a female-typical to a more male-typical pattern. 
We first identified CpG sites that were sexually di- 
morphic and defined them as those that displayed a sig- 
nificant difference in methylation levels between control 
females and control males (p value <0.05 measured by 
the Student's t test; FDR approximately 7% to 13%). This 
analysis identified about 12,000-20,000 sites in each 
brain region. For each site, we arbitrarily defined the 
male methylation level as 0 and the female level as 100. 
The XX + T methylation levels at these same sites were 
normalized to this scale and graphed on a continuum 
between 0 and 100. 

Almost all the sexually dimorphic CpG sites were 
more similar to females than males in the XX + T group 
at PN4 in both the striatum and the BNST/POA 
(Figure 4a,c); that is, testosterone treatment on the day 
of birth had not masculinized the methylome by PN4. 
However, by PN60, a large number of sexually dimorphic 
CpG sites in both the striatum and BNST/POA demon- 
strated methylation levels more similar to males than to 
females (p value <2.20E - 16, Kolmogorov-Smirnov test) 
(Figure 4b,d). The shift toward male values was more pro- 
nounced in the striatum than in the BNST/POA (compare 
Figure 4b to 4d). The distribution of sexually dimorphic 
CpG sites in PN60 XX + T resembles a bimodal distribu- 
tion (Figure 4b). These findings have two important impli- 
cations. The first is that neonatal testosterone treatment of 
females induces a clear shift to more masculine methylation 
patterns that is only apparent as the animal ages. The 
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3.43E- 


-03 


Hyperactive behavior 


3.47E- 


-03 


Schizophrenia 
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9.03E- 


-03 


Oligodendroglioma 


9.64E - 


-03 


Spina bifida 
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Table 3 Testosterone-driven methylation changes at several genes are maintained into adulthood 



Gene symbol 


A Me 
at PN4 


Fragment coordinate 
at PN4 


A Me 
at PN60 


Fragment coordinate 
at PN60 


Associated genomic features 
Prom GB CGI 


Striatum 
















4921515J06Rik 


0.22 


chr3:1 08742959-1 08743236 


-0.13 


chr3:1 08742959-1 08743236 


X 


X 


X 


Taf4b 


-0.11 


chrl 8:1 5048075-1 504831 8 


-0.1 


chrl 8:1 5048075-1 504831 8 




X 




F8a 


-0.1 


chrX:70473904-70474128 


-0.13 


chrX:70474305-70474419 




X 


X 


Fmr1 


-0.1 


chrX:65932427-65932537 


-0.11 


chrX:65932427-65932537 




X 


X 


Kcnql 


0.28 


chr7:1 50455734-1 50456023 


-0.14 


chr7:1 50481 372-1 50481 523 




X 


X 


Rbbp7 


-0.23 


chrX:1 591 98854-1 591 99070 


-0.12 


chrX:1 591 98688-1 591 98839 


X 


X 


X 


Sox3 


-0.11 


chrX:581 45676-581 45842 


-0.11 


chrX:581 46499-581 4661 9 


X 


X 


X 


Dab1 


-0.1 


chr4:1 04298501 -104298709 


-0.11 


chr4:1 04275571 -104275789 




X 




Nnat 


-0.13 


chr2:1 57385832-1 57386011 


-0.12 


chr2:1 57386045-1 5738621 4 


X 


X 


X 


Arid3b 


-0.19 


chr9:57685767-57685957 


-0.11 


chr9:57685767-57685957 


X 






Gripl 


-0.21 


chr10:1 19402892-1 19403155 


-0.12 


chr10:1 19402892-1 19403155 




X 


X 


Lonrf3 


-0.1 


chrX:33868422-33868652 


-0.12 


chrX:33869078-33869231 


X 


X 


X 


Clybl 


-0.15 


chrl 4:1 22662639-1 2266281 5 


-0.12 


chrl 4:1 22629995-1 226301 46 




X 




Sorcs2 


-0.13 


chr5:36720053-36720329 


-0.12 


chr5:365 1 1 861 -365 1 2058 




X 




2610018G03Rik 


-0.11 


chrX:481 94982-481 95 124 


-0.12 


chrX:481 94982-481 95 124 




X 


X 


Rap2c 


-0.17 


chrX:48370998-48371218 


-0.12 


chrX:48370998-48371218 


X 


X 


X 


Fndc3b 


0.12 


chr3:27382887-27383015 


-0.21 


chr3:2736871 6-27368995 




X 




Ubash3b 


-0.12 


chr9:40872064-40872268 


-0.14 


chr9:40872064-40872268 




X 




Foxkl 


-0.25 


chr5:142921051-142921281 


-0.27 


chr5:1 42921 051 -1 42921 281 




X 




BNST/POA 
















Igfbp7 


-0.13 


chr5:77809526-77809779 


-0.21 


chr5:77786342-77786589 




X 




Odz3 


-0.11 


chr8:49626485-49626755 


0.12 


chr8:4939501 2-49395267 




X 




Emd 


-0.13 


chrX:71 500275-71 500386 


-0.11 


chrX:71 500067-71 500242 


X 


X 


X 


Herc3 


-0.12 


chr6:58856760-58856872 


-0.12 


chr6:58856760-58856872 


X 


X 


X 


Commdl 


-0.16 


chrl 1:22873668-22873935 


-0.12 


chrl 1:22872579-22872758 


X 


X 


X 


Bcor 


-0.14 


chrX:1 1715730-1 1715985 


-0.11 


chrX:1 1703662-1 1703859 




X 




Napl 15 


-0.12 


chr6:58856760-58856872 


-0.12 


chr6:58856760-58856872 


X 


X 


X 


Gpr179 


-0.14 


chrl 1:971 93837-971 941 08 


-0.11 


chrl 1:971 97795-971 97992 




X 




Zrsrl 


-0.16 


chrl 1:22873668-22873935 


-0.12 


chrl 1:22872579-22872758 


X 


X 


X 


Lonrf3 


-0.12 


chrX:33868422-33868652 


-0.1 


chrX:33869078-33869231 


X 


X 


X 


Sdk1 


-0.13 


chr5:142590169-142590448 


0.12 


chr5:1 423 1 241 0-1 4231 2648 




X 





Differential CpG methylation values between XX and XX + T at PN4 and PN60 are shown. Negative delta methylation indicates higher methylation in XX + T. 
Genomic features associated with the differentially methylated region are indicated by the last three columns. An 'X' indicates that the fragment overlaps with a 
feature {Prom promoter, GB gene body, CGI CpG island). 



second is that the methylome of the striatum may be more 
responsive to the effects of neonatal testosterone than the 
BNST/POA. Additional file 7: Table S6 contains the full list 
of masculinized genes. 

Expression analysis 

We also carried out transcriptome profiling on five bio- 
logical replicates per experimental group using cDNA 
microarrays. Each replicate consisted of a unique pool of 
tissue from three animals from the same set of samples 



used for methylomic analysis. After filtering for high- 
quality expression data, 57 samples and 13,776 array 
probes were retained for differential expression analysis. 
According to principal component analysis (PCA), age and 
brain region explained 78% and 48% of the expression vari- 
ance, respectively, while experimental group differences 
only explained up to 0.05% (Figure 5a; Additional file 8: 
Table S7). 

To identify genes with sex-specific expression pat- 
terns and those in which expression is regulated by 
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Methylation Scale 

Figure 4 CpG methylation patterns are more masculine in XX + T than XX at PN60. Mean XX + T methylation of the sites that display sizable basal 
sex differences (delta methylation (X x-xY) ^15%, p value <0.05) are plotted on a scale where 0 corresponds to XY and 100 corresponds to XX methylation 
levels. At the bottom of each histogram the first number in the box represents the number of masculinized sites (methylation scale <50) divided by the 
second number which is the number of total sites in the histogram. PN4 (a) and PN60 (c) striatum, PN4 (b), and PN60 BNST/POA (d). In both brain regions, 
the PN60 distribution was significantly different from PN4 ip value <2.20E- 16, Kolmogorov-Smirnov test). 



testosterone, we tested each of the possible pairwise 
comparisons (XX vs. XY, XX vs. XX + T, XX + T vs. 
XY) in different ages and brain regions. Significant 
differential expression signals were only observed for 
(1) XX vs. XY comparison at PN4 in BNST/POA (n = 5 
per group, Benjamini-Hochberg (BH)-adjusted p value 
(FDR) <0.05, 21 genes) and (2) XX vs. XX + T at PN60 in 
the striatum (n = 5 per group, BH-adjusted p value (FDR) 
<0.05, 99 genes) (Figure 5b,c). qRT-PCR was performed to 
validate the microarray results. All genes detected as sig- 
nificantly different in the BNST/POA XX vs. XY PN4 
microarray showed significant differences by PCR in the 
expected direction (Additional file 5: Figure S2). We did 
not have the biological material to validate the striatal XX 
vs. XX + T PN60 microarray. 

None of the genes in the striatal PN4 XX vs. XX + T 
comparison survived BH adjustment. However, by day 
60 of life, 99 genes demonstrated statistically significant 
testosterone-dependent expression changes (BH-adjus- 
ted p value (FDR) <0.05). This finding was consistent 
with the trend observed in the methylation data. In 
addition, the gene expression differences associated with 
testosterone treatment were less pronounced in the BNST/ 



POA than in the striatum in adulthood (Additional file 8: 
Table S7). IPA analysis of differentially expressed genes in 
XX vs. XX + T PN60 striatum showed enrichment of gene 
categories such as formation ofactin filaments, formation of 
neurites, myelination of neurites, long-term potentiation, 
and abnormal morphology of dopaminergic neurons 
(Additional file 9: Table S8). 

We next tested whether differentially expressed genes 
also showed differential methylation. We identified a 
number of genes where the CpG methylation state and 
gene expression patterns correlated in the XX vs. XX + T 
comparison in PN60 striatum. Of the 99 striatal genes 
that were differentially expressed between XX and XX + T 
at PN60, 10 genes also showed differential methylation 
(Table 4). There were no genes that showed both differen- 
tial expression and methylation in the BNST/POA. 

Overall, we identified genes and pathways that are sub- 
ject to regulation by testosterone at the mRNA level. Our 
data suggest that testosterone regulates gene transcription 
in a highly context-dependent manner (e.g., the effect of 
testosterone appeared to be stronger in the striatum). 
Finally, although some level of correlation was ob- 
served between methylation and expression changes, 



Ghahramani et al. Biology of Sex Differences 2014, 5:8 
http://www.bsd-journal.eom/content/5/1/8 



Page 14 of 18 



a 



Cluster dendrogram with AU/BP values (%) 



Age 



Distance: euclidean 
Cluster method: average 



ic\nn co co c\j n * \- >in w X 1/1 T . 




C.JJ o«H 
coZX 



CO 5 '.£= i= - i; o C/3 C/) (/)(- 
"- </>«(/)(/) CO I ItOHI- 
I iZ0007mm 

CLCL 



X t XX< 5vv Szxx 

l¥ XXO XX , X | ZCOX 

_ kt< !Q -<« CL 0<0 

"■^OOl-OOO <oocl 

mw ZZ oZzz 2m 

-'££<ocococo £ 



Q-(D^^ ZZZ 
"^qIqI Q-Q-Q-j 



Jl_ 



Brain region I 



b P value distributions C Top 20 differentially expressed genes 



PN4 BNST/POA XX vs. XY 




-1.2 -0.8 -0.4 0.0 0.4 0.8 1.2 




I I I I I I I I I I I I I I I 
-1.2 -0.8 -0.4 0.0 0.4 0.8 1.2 



Figure 5 RNA samples cluster based on age and brain region differences, (a) Dendrogram representing hierarchical clustering of samples 
based on 1 -Pearson correlation of the genome-wide expression profiles. Red values are approximately unbiased {AU) p values, and green values are 
bootstrap probability {BP) values. Values on the edges of the dendrogram are p values (%). Samples tended to cluster by age (PN4/PN60) and 
brain region (striatum/BNST/POA). (b) p value distributions of comparison groups that demonstrated a significant differential expression signal 
(PN4 XX vs. XY BNST/POA and PN60 XX vs. XX + T striatum), (c) Each bar plot represents the log 2 fold changes in gene expression of the top 20 
upregulated and downregulated probe sets. For each plot the relative gene expression fold changes are shown for the XY or XX + T relative to 
the XX. Positive fold change values represent increased expression in XY or XX + T relative to XX mice. 



the genes identified by these two analyses were largely 
non-overlapping. The ability of methylation differ- 
ences to predict expression changes at the single gene 
level was limited. 



Discussion 

Our study shows that the organizational effect of testos- 
terone on the brain is late-emerging and markedly modi- 
fies the epigenetic CpG methylation landscape of the 
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Table 4 Genes that show testosterone-driven changes in 
both CpG methylation and gene expression 



Gene 


L092 fold 
change 


Expression difference 
p value 


i\ ivie at 
PN60 


Ppp2r5c 


1.07 


9.66E - 05 


0.10 


Nrxnl 


0.95 


3.63E-04 


0.10 


Lpinl 


0.71 


8.18E-05 


0.10 


Eglnl 


0.53 


2.41 E- 04 


0.13 


Binl 


0.47 


2.08E - 04 


0.11 


Inpp5a 


0.42 


3.80E - 04 


-0.13 


Fhll 


0.40 


4.1 1 E - 05 


0.15 


Mras 


-0.42 


3.09E - 04 


0.13 


E130203B14Rik 


-0.61 


2.75E-04 


0.12 


Wftid 


-0.66 


1 .26E - 04 


0.10 



Genes that show testosterone-driven changes in both CpG methylation and 
gene expression in the PN60 XX vs. XX + T comparison in the striatum. Nega- 
tive values in the second column indicate lower expression in XX + T relative to 
XX. In the fourth column, negative values indicate lower methylation in XX + T 
relative to XX. 



brain. We have established the first genome-wide and 
quantitative map of testosterone-induced CpG methyla- 
tion changes in two sexually dimorphic brain regions. 
We show that the majority of methylation changes 
caused by testosterone are not established in the first 
few days following the initial exposure to testosterone. 
Instead, testosterone's effects are most obvious in adult- 
hood where it has masculinized a large number of 
sexually dimorphic CpG sites, which show increased 
methylation in XX + T and XY compared to XX. These 
results run counter to our initial prediction that sex dif- 
ferences in methylation are established during early de- 
velopment, most of which remain stable over time. Since 
the organizational effects of testosterone appear much 
later in life, organization by testosterone may occur via 
early programming on relatively few genes. This small 
initial effect is what sets up the brain to respond in a 
particular fashion to other events during postnatal 
development. Since all animals in the study were gonad- 
ectomized before puberty and provided with identical 
hormone replacement prior to sacrifice, differences seen 
during adulthood likely reflect programming by perinatal 
testosterone. However, the full manifestation of these 
differences may require the presence of testosterone in 
adulthood. 

All animals were given T implants at PN45 to cap- 
ture the full effect of perinatal organization. Previous 
studies had shown that the full manifestation of some 
organizational effects requires activational hormonal effects 
(in other words, adult levels of circulating hormone levels). 
For instance, female rats given testosterone perinatally will 
not show increased mounting behavior unless they are also 
given a dose of testosterone in adulthood [38]. Therefore, 



we hypothesized that the effects of perinatal T on the 
methylome would not fully manifest themselves unless 
adult circulating hormonal levels were also present. If 
no hormone replacement had taken place, it is likely 
that we would have found fewer differences between 
XX and XX + T. It is also likely that the XX + T ani- 
mals would have experienced less masculinization of 
their methylomes. 

Another unexpected finding with respect to the testos- 
terone effects is that the shift toward a male-like pattern 
of CpG methylation in XX + T mice was more pro- 
nounced in the striatum than in the BNST/POA. This 
was surprising since the BNST/POA displays some of the 
most prominent anatomical and neurochemical sex 
differences in the brain that result from organization 
by gonadal hormones, whereas sex differences in the 
striatum are comparatively modest in terms of neuro- 
anatomy. Our data therefore show that neuroanatomical 
and neurochemical markers of sex differences may not 
fully reflect the sensitivity of a brain region to gonadal 
hormones. 

We also identified sets of testosterone-regulated loci 
that clearly maintained differences in 5-mC from PN4 to 
PN60 in both the striatum and BNST/POA, although 
these were a small minority. The overwhelming majority 
of testosterone-affected loci showed dynamic CpG me- 
thylation patterns. While this is not in agreement with 
the classic view of DNA methylation as a permanent 
epigenetic mark, our data are consistent with the find- 
ings of several recent studies demonstrating develop- 
mental or experience-dependent regulation of 5-mC in 
the mammalian brain [39-43]. For example, Schwarz 
et al. have shown that sex differences in methylation pat- 
terns at the Esrl, Esr2, and Pr promoters are dynamic 
across the life span [15]. Although, we did not identify 
the methylation patterns of Esrl, Esr2> and Pr promoters 
as significantly influenced by sex or testosterone expos- 
ure, we were not surprised given the different species 
used in each study (mice in our study and rats in 
Schwarz et al), the difference in tissue examined (both 
the BNST and POA here vs. just the POA in Schwarz 
et al), and the fact that we required our differentially 
methylated fragments to show consistent methylation 
changes in several adjacent CpG sites (Schwarz et al. 
considered each CpG site on its own). 

When we examined the characteristics of the genes as- 
sociated with the testosterone-modified CpGs, we found 
a significant enrichment of genes that are expressed in 
the brain, particularly those involved in synaptic function. 
This suggests that the effect of the neonatal testosterone is 
not random and that testosterone specifically alters 
the methylation of neural-related genes. Furthermore, 
the effects of testosterone (particularly on X-linked genes) 
appear to be brain region-specific. Interestingly, among 
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the genes that were differentially methylated in the stri- 
atum due to testosterone, a substantial number encoded 
signaling components associated with increased or de- 
creased risk of Parkinson's disease. We found that these 
genes overlapped with genes that show strong sex differ- 
ences in expression in dopaminergic neurons from Parkin- 
son's disease patients [44]. 

In the BNST/POA, sexually differentiated rates of 
apoptosis (male < female) driven by testosterone expos- 
ure is one of the major events leading to the volumetric 
sex differences in this region during the sensitive period 
[22,45]. Consistent with this, we found genes involved 
in apoptosis in the testosterone-affected dataset at 
both PN4 and PN60. Perhaps most intriguing, we had 
genes related to the proliferation of neuronal cells in 
our dataset. Recent evidence has shown that the 
maintenance of sexual dimorphism resulting from the 
organizational effects of testosterone may require rein- 
forcement in the form of pubertal hormones [36,46,47]. 
At least some of this reinforcement may take the 
form of sexually differentiated cell addition in several 
sexually dimorphic brain regions including the POA 
[36,48,49]. 

In studying the genome-wide methylation profiles of 
the striatum and BNST/POA both neonatally and later 
in life, we found that the methylation patterns of a large 
number of genes differed between the two sexes. Male- 
to-female comparisons displayed a marked enrichment 
of methylation on the X chromosome in females which 
could potentially be accounted for by X chromosome in- 
activation. Most of the X-linked genes that showed sex- 
specific differences in CpG methylation were found in 
both brain regions. This strongly suggests that the 
methylation differences we are detecting are real as the 
process of X inactivation affects X-linked genes similarly 
across tissues [34]. We also identified widespread auto- 
somal gene methylation differences between males and 
females. In contrast to the results from X-linked genes, 
the vast majority of autosomal genes that show sex dif- 
ferences were unique to either the striatum or the 
BNST/POA. Most of these genes were more methylated 
in males. One possible explanation is that the inactivated 
X chromosome could potentially act as a 'heterochroma- 
tin sink,' sequestering the factors required for gene silen- 
cing (DNA methyltransferases, etc.) at autosomal loci 
[50-52]. 

Most interestingly, concordant with the findings from 
our analysis of CpG methylation in the striatum, we 
found that testosterone's effect on gene expression was 
late-emerging. Finally, when we explored the extent to 
which changes in methylation levels contributed to gene 
expression differences, we found evidence for correla- 
tions between gene expression patterns and methylation 
profiles at some genes associated with DMRs, but in 



general, methylation differences did not predict differ- 
ences in gene expression. This observation may not be 
surprising given that we worked with tissues that were 
relatively heterogeneous. In addition, many CpG methy- 
lation differences may not be associated with gene ex- 
pression changes. Lyko et al. showed that there is a 
strong correlation between CpG methylation and spli- 
cing sites including those that have the potential to yield 
alternative exons [53]. Therefore, an interesting avenue 
of research will be determining whether testosterone- 
induced methylation differences can specifically regulate 
splicing, rather than transcription. 

There are several limitations to this study. This 
work represents a snapshot of the DNA methylation 
landscape at two ages, while the brain may display a 
vast array of epigenetic states as it passes through different 
stages of development. Longitudinal study designs exam- 
ining DNA methylation at different life stages could pro- 
vide a comprehensive picture of how the epigenome 
is modified over time. Because testosterone is known 
to alter the proportion of specific cell types comprising 
the BNST/POA and possibly the striatum, group differ- 
ences in levels of methylation reported here could be the 
result of testosterone- or sex-specific regulation of the cell 
types in the dissected tissue, which in turn differ in 
their methylation of specific genes, or could reflect 
direct testosterone effects on the methylome of cell 
types common to the different groups. In addition, 
DNA methylation is associated with other epigenetic 
alterations, especially histone modifications and RNAi 
pathways. Studies of these other epigenetic changes 
are crucial to identifying common mechanisms under- 
lying sex differences in epigenetic regulation. Finally, 
different brain regions are expected to display differ- 
ent epigenetic marks across their genomes, and epi- 
genetic profiling across functionally discrete brain 
areas will be important in identifying tissue-specific sex 
differences. 



Conclusions 

Taken together, our results suggest that early testoster- 
one exposure has broad effects on brain methylation 
patterns particularly during adulthood and that the 
emergence of sex differences in the brain may be a grad- 
ual process that is cemented over the organism's life. 
Our data provide a new perspective by showing that 
most sex differences in CpG methylation are dynamic 
and not the result of acute modifications in response to 
hormones. Clearly, additional studies of genome-scale 
methylation maps will be important to give us a full 
understanding of the long-lasting influences of early 
hormone exposure on DNA methylation dynamics of 
the brain. 
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Additional files 



Additional file 1: Table SI. Summary information for all individual 
striatum and BNST/POA samples included in this study. To enable our 
comparison, RRBS libraries were sequenced and aligned to the mouse 
mm9 genome. Twelve striatum RRBS and twelve BNST/POA RRBS libraries 
(125 million total reads, average coverage of 58x, and 46% mapability) 
were chosen for comparison, and summary statistics are shown. 
Genome coverage represents the number of CpG sites in the genome 
calculated on both strands. Fragment coverage reflects the number 
of CpG sites in the fragments of sizes 1 00 to 300. Data coverage 
is the number of CpG sites covered by at least one read, and 
analysis coverage refers to the number of CpG sites covered by 
at least four reads. 

Additional file 2: Table S2. Genes demonstrating differential 
methylation in XX vs. XY animals. Differential methylation is defined by a 
z-score threshold that results in a false discovery rate less than 10% and 
DNA methylation difference >10% as represented by delta methylation. 
Negative delta methylation values represent increased methylation 
in XX. 

Additional file 3: Table S3. Differentially methylated genes in the 
testosterone dataset (XX vs. XX + T). Differential methylation is defined by 
a z-score threshold that results in a false discovery rate less than 10% and 
DNA methylation difference >10% as represented by delta methylation. 
Negative delta methylation values represent increased methylation levels 
in XX + T relative to XX mice. 

Additional file 4: Table S4. Genes whose methylation patterns were 
different between the XX and XX + T and between the XX and XY 
animals. Differential methylation is defined by a z-score threshold that 
results in a false discovery rate less than 10% and DNA methylation 
difference >10% as represented by delta methylation. Negative delta 
methylation values in the XX vs. XX + T columns represent increased 
methylation in XX + T relative to XX. In the XX vs. XY dataset, negative 
delta methylation values represent increased methylation in XX 
relative to XY. 

Additional file 5: Validation of RRBS and microarray results. 
Figure SI. Validation of RRBS data by traditional (Sanger) bisulfite 
sequencing. One locus from each brain region that showed both 
testosterone-affected and sex-specific methylation differences at PN60 was 
selected for validation (for the striatum, Micalll; BNST/POA: Fzd9). n = 2-3 per 
group. Figure S2. qPCR validation of genes that were detected as significantly 
differentially expressed between XX and XY in the BNST/POA at PN4. 
Error bars represent the standard error of the mean from 3 to 4 
biological replicates from each group. Expression is relative to GAPDH 
and is normalized to XX. 

Additional file 6: Table S5. Ingenuity Pathway Analysis (IPA, www. 
ingenuity.com) was performed on the entire dataset of genes whose 
methylation was affected by testosterone at PN60 in striatum and BNST/ 
POA. The IPA networks, functions, and canonical pathways enriched in 
this dataset are provided in separate tabs. 

Additional file 7: Table S6. Genes associated with CpG sites where 
methylation was masculinized by testosterone in the XX + T group. 
Genes were deemed masculinized if they contained at least one CpG site 
in either the gene body or the promoter demonstrating more male 
typical methylation patterns. 

Additional file 8: Table S7. Lists of differentially expressed genes 
between XX vs. XX + T and between XX vs. XY in striatum and BNST/POA 
at both time points (n = 5 per group, Benjamini-Hochberg-adjusted p 
value <0.05). Negative fold change values represent decreased expression 
in XX + T or XY relative to XX mice. 

Additional file 9: Table S8. Ingenuity Pathway Analysis of genes 
associated with XX vs. XX + T comparison at PN60 in the striatum. The 
IPA networks, functions, and canonical pathways enriched in this dataset 
are provided in separate tabs. 
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